The purpose of these benchmarks is to compare several integrators for use in molecular dynamics simulation. We will use a simulation of liquid argon form the examples of NBodySimulator as test case.

using ProgressLogging
using NBodySimulator, OrdinaryDiffEq, StaticArrays
using Plots, DataFrames, StatsPlots

function setup(t)
    T = 120.0 # K
    kb = 1.38e-23 # J/K
    ϵ = T * kb # J
    σ = 3.4e-10 # m
    ρ = 1374 # kg/m^3
    m = 39.95 * 1.6747 * 1e-27 # kg
    N = 216
    L = (m*N/ρ)^(1/3)
    R = 2.25σ
    v_dev = sqrt(kb * T / m) # m/s

    _L = L / σ
    _σ = 1.0
    _ϵ = 1.0
    _m = 1.0
    _v = v_dev / sqrt(ϵ / m)
    _R = R / σ

    bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
    lj_parameters = LennardJonesParameters(_ϵ, _σ, _R)
    pbc = CubicPeriodicBoundaryConditions(_L)
    lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters));
    simulation = NBodySimulation(lj_system, (0.0, t), pbc, _ϵ/T)

    return simulation
end

setup (generic function with 1 method)

In order to compare different integrating methods we will consider a fixed simulation time and change the timestep (or tolerances in the case of adaptive methods).

function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt, b, gc, memalloc = @timed solve(prob, alg();
            save_everystep=false, progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = ΔE
        rts[alg] = rt
        bytes[alg] = b
        allocs[alg] = memalloc
        nt[alg] = sol.destats.naccept
        nf[alg] = sol.destats.nf + sol.destats.nf2
    end
end

function run_benchmark!(results, t, integrators, tol...; c=ones(length(integrators)))
    @progress "Benchmark at t=$t" for τ in zip(tol...)
        runtime = Dict()
        ΔE = Dict()
        nt = Dict()
        nf = Dict()
        b = Dict()
        allocs = Dict()
        cfg = config(integrators, c, τ...)

        GC.gc()
        benchmark(ΔE, runtime, b, allocs, nt, nf, t, cfg)
        get_tol(idx) = haskey(cfg[idx], :dt) ? cfg[idx].dt : (cfg[idx].abstol, cfg[idx].rtol)

        for (idx,i) in enumerate(integrators)
            push!(results, [string(i), runtime[i], get_tol(idx)..., abs(ΔE[i]), nt[i], nf[i], c[idx]])
        end
    end
    return results
end

run_benchmark! (generic function with 1 method)

We will consider symplectic integrators first

symplectic_integrators = [
    VelocityVerlet,
    # VerletLeapfrog,
    PseudoVerletLeapfrog,
    McAte2,
    # CalvoSanz4,
    # McAte5,
    Yoshida6,
    KahanLi8,
    # SofSpa10
]

5-element Array{DataType,1}:
 VelocityVerlet
 PseudoVerletLeapfrog
 McAte2
 Yoshida6
 KahanLi8

Let us run a short simulation to see the cost per timestep for each method

config(integrators, c, τ) = [ (alg=a, dt=τ*cₐ) for (a,cₐ) in zip(integrators, c)]

t = 35.0
τs = 1e-3

# warmup
c_symplectic = ones(length(symplectic_integrators))
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(symplectic_integrators, c_symplectic, τs))

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, symplectic_integrators, τs)

5 rows × 7 columns

	integrator	runtime	τ	EnergyError	timesteps	f_evals	cost
	String	Float64	Float64	Float64	Int64	Int64	Float64
1	VelocityVerlet	58.4098	0.001	0.00366562	35000	70002	1.0
2	PseudoVerletLeapfrog	59.2155	0.001	0.0208884	35000	105002	1.0
3	McAte2	57.6713	0.001	0.00980374	35000	105002	1.0
4	Yoshida6	229.896	0.001	0.0219439	35000	525002	1.0
5	KahanLi8	530.512	0.001	0.0223923	35000	1225002	1.0

The cost of a timestep can be computed as follows

c_symplectic .= results[!, :runtime] ./ results[!, :timesteps]
c_Verlet = c_symplectic[1]
c_symplectic /= c_Verlet

5-element Array{Float64,1}:
 1.0
 1.0137927860016844
 0.9873557593070432
 3.935903166715239
 9.082574686069826

were we have normalized the cost to the cost of a VelocityVerlet step.

Let us now benchmark the solvers for a fixed simulation time and variable timestep

t = 40.0
τs = 10 .^range(-4, -3, length=5)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)

25 rows × 7 columns

	integrator	runtime	τ	EnergyError	timesteps	f_evals	cost
	String	Float64	Float64	Float64	Int64	Int64	Float64
1	VelocityVerlet	679.942	0.0001	0.00279731	400000	800002	1.0
2	PseudoVerletLeapfrog	671.562	0.000101379	0.000301274	394558	1183676	1.01379
3	McAte2	683.144	9.87356e-5	0.000509045	405123	1215371	0.987356
4	Yoshida6	684.565	0.00039359	0.0128874	101629	1524437	3.9359
5	KahanLi8	665.012	0.000908257	0.0357851	44041	1541437	9.08257
6	VelocityVerlet	378.816	0.000177828	0.00394572	224937	449876	1.0
7	PseudoVerletLeapfrog	375.186	0.000180281	0.00147033	221877	665633	1.01379
8	McAte2	386.271	0.000175579	0.00185951	227818	683456	0.987356
9	Yoshida6	382.206	0.000699914	0.0289537	57150	857252	3.9359
10	KahanLi8	381.662	0.00161514	0.00651373	24766	866812	9.08257
11	VelocityVerlet	212.669	0.000316228	0.0134205	126492	252986	1.0
12	PseudoVerletLeapfrog	214.316	0.000320589	0.0113341	124771	374315	1.01379
13	McAte2	217.365	0.000312229	0.00327997	128111	384335	0.987356
14	Yoshida6	218.757	0.00124464	0.013221	32138	482072	3.9359
15	KahanLi8	209.768	0.00287216	0.0680963	13927	487447	9.08257
16	VelocityVerlet	121.405	0.000562341	0.00608374	71132	142266	1.0
17	PseudoVerletLeapfrog	119.33	0.000570098	0.00520789	70164	210494	1.01379
18	McAte2	124.652	0.000555231	0.00932556	72043	216131	0.987356
19	Yoshida6	121.941	0.00221332	0.000567354	18073	271097	3.9359
20	KahanLi8	119.632	0.00510751	0.181431	7832	274122	9.08257
21	VelocityVerlet	67.8505	0.001	0.0110783	40001	80004	1.0
22	PseudoVerletLeapfrog	66.9239	0.00101379	0.0224749	39456	118370	1.01379
23	McAte2	68.7452	0.000987356	0.00814607	40513	121541	0.987356
24	Yoshida6	69.461	0.0039359	0.0120663	10163	152447	3.9359
25	KahanLi8	67.6691	0.00908257	0.308016	4405	154177	9.08257

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

Looking at the runtime as a function of timesteps, we can observe that we have a linear dependency for each method, and the slope is the previously computed cost per step.

@df results plot(:timesteps, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Number of timesteps", ylabel="Runtime (s)")

We can also consider a longer simulation time

t = 100.0

τs = 10 .^range(-4, -3, length=5)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
#run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)

The energy error as a function of runtime is given by

#@df results plot(:EnergyError, :runtime, group=:integrator,
#    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

We can also look at the energy error history

function benchmark(energyerr, rts, ts, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE(t) = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = [ΔE(t) for t in sol.t[2:end]]
        rts[alg] = rt
        ts[alg] = sol.t[2:end]
    end
end

ΔE = Dict()
rt = Dict()
ts = Dict()
configs = config(symplectic_integrators, c_symplectic, 2.3e-4)
benchmark(ΔE, rt, ts, 45., configs)

plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:topleft);
for c in configs
    plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s", yscale=:log10)
end
plt

Now, let us compare some adaptive methods

adaptive_integrators=[
    # DPRKN
    DPRKN6,
    # DPRKN8,
    # DPRKN12,
    # others
    Tsit5,
    # Vern7,
    # Vern9
]

config(integrators, c, at, rt) = [ (alg=a, abstol=at, rtol=rt) for a in integrators]

t = 35.0
ats = 10 .^range(-20, -14, length=10)
rts = 10 .^range(-20, -14, length=10)

# warmup
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(adaptive_integrators, 1, ats[1], rts[1]))

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, adaptive_integrators, ats, rts)

20 rows × 8 columns

	integrator	runtime	abstol	reltol	EnergyError	timesteps	f_evals	cost
	String	Float64	Float64	Float64	Float64	Int64	Int64	Float64
1	DPRKN6	32.9054	1.0e-20	1.0e-20	0.0363926	6035	45238	1.0
2	Tsit5	18.1985	1.0e-20	1.0e-20	0.1444	3281	21717	1.0
3	DPRKN6	32.995	4.64159e-20	4.64159e-20	0.00214219	6064	45693	1.0
4	Tsit5	18.4353	4.64159e-20	4.64159e-20	0.287378	3262	21687	1.0
5	DPRKN6	31.6542	2.15443e-19	2.15443e-19	0.0618527	5916	44237	1.0
6	Tsit5	18.8475	2.15443e-19	2.15443e-19	0.942605	3317	22281	1.0
7	DPRKN6	31.9917	1.0e-18	1.0e-18	0.012243	5948	44573	1.0
8	Tsit5	17.959	1.0e-18	1.0e-18	0.798038	3231	21477	1.0
9	DPRKN6	32.1473	4.64159e-18	4.64159e-18	0.0302664	5992	44951	1.0
10	Tsit5	18.4408	4.64159e-18	4.64159e-18	0.995148	3236	21531	1.0
11	DPRKN6	31.5106	2.15443e-17	2.15443e-17	0.056756	5894	44062	1.0
12	Tsit5	18.504	2.15443e-17	2.15443e-17	0.732307	3281	21957	1.0
13	DPRKN6	30.4919	1.0e-16	1.0e-16	0.0582733	5767	42690	1.0
14	Tsit5	18.0888	1.0e-16	1.0e-16	0.805652	3270	21687	1.0
15	DPRKN6	32.9233	4.64159e-16	4.64159e-16	0.0369794	5928	44475	1.0
16	Tsit5	18.2209	4.64159e-16	4.64159e-16	0.246302	3256	21621	1.0
17	DPRKN6	31.928	2.15443e-15	2.15443e-15	0.0209199	5917	44076	1.0
18	Tsit5	18.6343	2.15443e-15	2.15443e-15	0.0214593	3269	21729	1.0
19	DPRKN6	31.4097	1.0e-14	1.0e-14	0.100083	5872	43768	1.0
20	Tsit5	18.3741	1.0e-14	1.0e-14	0.314547	3259	21783	1.0

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

If we consider the number of function evaluations instead, we obtain

@df results plot(:EnergyError, :f_evals, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Number of f evals")

We can also consider a longer simulation time

t = 100.0

ats = 10 .^range(-20, -14, length=10)
rts = 10 .^range(-20, -14, length=10)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
#run_benchmark!(results, t, integrators, ats, rts)

The energy error as a function of runtime is given by

#@df results plot(:EnergyError, :runtime, group=:integrator,
#    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

---